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ABSTRACT 

We describe a method for the numerical computation of the propagation of primary and secondary nucleons, 
primary electrons, and secondary positrons and electrons. Fragmentation and energy losses are computed using 
realistic distributions for the interstellar gas and radiation fields, and diffusive reacceleration is also incorporated. 
The models are adjusted to agree with the observed cosmic -ray B/C and 10 Be/ 9 Be ratios. Models with diffusion 
and convection do not account well for the observed energy dependence of B/C, while models with reacceleration 
reproduce this easily. The height of the halo propagation region is determined, using recent 10 Be/ 9 Be measure- 
ments, as > 4 kpc for diffusion/convection models and 4-12 kpc for reacceleration models. For convection 
models we set an upper limit on the velocity gradient of dV /dz < 7 km s" 1 kpc" 1 . The radial distribution of 
cosmic -ray sources required is broader than current estimates of the SNR distribution for all halo sizes. Full 
details of the numerical method used to solve the cosmic -ray propagation equation are given. 

Subject headings: cosmic rays — diffusion — elementary particles — Galaxy: general — ISM: abundances — 
ISM: general 



1. INTRODUCTION 

A numerical method and corresponding computer code for 
the calculation of Galactic cosmic-ray propagation has been de- 
veloped, which is a further development of the approach de- 
scribed by Strong & Youssefi (1995) and Strong (1996). Pri- 
mary and secondary nucleons, primary and secondary elec- 
trons, and secondary positrons are included. The basic spa- 
tial propagation mechanisms are (momentum-dependent) dif- 
fusion and convection, while in momentum space energy loss 
and diffusive reacceleration are treated. Fragmentation and en- 
ergy losses are computed using realistic distributions for the 
interstellar gas and radiation fields. Preliminary resu lts were 
presented in Strong & Moskalenko (1997) (hereafter Paper I) 



and full results for protons, Helium, posit rons, and electrons 
i n Moskalenko & Strong (1998a) (hereafter |Paper Il|). In P aper 
II we referred the description of the numerical method to the 
present paper (Paper III), and full details are now given. Re- 
sults for gamma-rays and synchrotron rad iation wil l be given in 



Moskalenko & Strong (1998b) (hereafter ^aper IV|) 



pei 

We note that our positron predictions from [Paper II| have been 



compared with more recent absolute measurements in Barwick 
et al. (1998) and the agreement is good; for the positrons this 
new comparison has the advantage of being independent of the 
electron spectrum, unlike t he positron/electron ratio which was 
the main focus of paper 5 The ultimate goal is to combine all 
constraints includi ng gamm a-ray and synchrotron spectra; this 
will be pursued in paper IV . 

The rationale for our approach was given previously (P aper 
I, paper lip . Briefly, the idea is to develop a model which simul- 
taneously reproduces observational data of many kinds related 
to cosmic-ray origin and propagation: directly via measure- 
ments of nuclei, electrons, and positrons, indirectly via gamma 
rays and synchrotron radiation. These data provide many in- 
dependent constraints on any model and our approach is able 
to take advantage of this since it must be consistent with all 
types of observation. We emphasize also the use of realistic 
astrophysical input (e.g. for the gas distribution) as well as the- 



oretical developments (e.g. reacceleration). The code is suf- 
ficiently general that new physical effects can be introduced as 
required. We aim for a 'standard model' which can be improved 
with new astrophysical input and additional observational con- 
straints. For interested users our model is available in the public 
domain on the World Wide Web. 

It was poi nted out many years ago (see G inzburg, Khazan & 
Ptuskin 1980, |Berezinskii et al. 199C| ) that the interpretation of 
radioactive cosmic-ray nuclei is model-dependent and in partic- 
ular that halo models lead to a quite different physical picture 
from homogeneous models. The latter show simply a rather 
lower average matter density than the l ocal Galactic hydroge n 
(e.g., |Simpson & Garcia-Munoz 1988 , Lukasiak et al. 1994a), 
but do not lead to a meaningful estimate of the size of the con- 
finement region, and the correponding cosmic -ray 'lifetime' is 
model-dependent. In such treatments the lifetime is combined 
with the grammage to yield an 'average density'. For exam- 
ple Lukasiak et al. (1994a) find an 'average density' of 0.28 
cm" 3 compared to the local interstellar value of about 1 cm" 3 , 
indicating a z-extent of less than 1 kpc compared to the several 
kpc found in diffusive halo models. In the present work we use 
a model which includes spatial dimensions as a basic element, 
and so these issues are automatically addressed. 

The possible role of convection was shown by Jokipii 
(1976), and Jones (1979) pointed out its effect on the energy- 
dependence of the secondary /primary ratio. Recent papers give 
estimates for the halo size and limits on convection based on 



existing calculations (e.g., Webber, Lee & Gupta 1992), and in 



the present work we attempt to improve on these models with a 
more detailed treatment. 

Previous approaches to the spatial nucleon propagation prob- 
lem have been mainly analytical: Jones (1979), Freedman et 
al. (1980), Berezinskii et al. (1990), Webber, Lee & Gupta 
(1992) and Bloemen et al. (1993) treated diffusion/convection 
models in this way. A problem here is that energy losses are 
difficult to treat and in fact were apparently not included ex- 
cept by Webber, Lee & Gupta (1992), however even there not 
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explicitly. Bloemenetal. (1993) used the 'grammage' formula- 
tion rather than the explicit isotope ratios, and their propagation 
equation implicitly assumes identical distributions of primary 
and secondary source functions. These papers did not attempt 
to fit the low-energy (< 1 GeV/nucleon) B/C data (which we 
will show leads to problems) and also did not consider reac- 
celeration. It is clear than an analytical treatment quickly be- 
comes limited as soon as more realistic models are desired, and 
this is the main justification for the numerical approach pre- 
sented in this paper. The case of electrons and positrons is even 
more intracta ble analytically, although fairly general cases have 
been treated ( Lerche & Schlickeiser 1982[ ). Owens & Jokipii 
(1977a,b) adopted an alternative approach with Monte-Carlo 
simulations, for both nucleons and electrons. Recently Porter & 
Protheroe (1997) made use of this method for electrons. Both 
these applications are for 1-D propagation, in the z-direction 
only. This method allows realistic models to be computed, 
but would be very time-consuming for 2- or 3-D cases. Our 
method, using numerical solution of the propagation equation, 
is a practical alternative. Since most of these studies were done, 
the data on both stable and radioactive nuclei has improved con- 
siderably and thus a re-evaluation is warranted. 

Reaccelera tion has previously been handled using leaky -box 
calc ulations (|Letaw, Silberberg & Tsao 1993L Seo & Ptuskin 
1994, fteinbach & S imon 1995[ ); this has the advantage of al- 
lowing a full reaction network to be used (far beyond what is 
possible in the present approach), but suffers from the usual 
limitations of leaky -box models, especially concerning radioac- 
tive nuclei, which were not included in these treatments. Our 
simplified reaction network is necessary because of the added 
spatial dimensions, but we believe it is fully sufficient for our 
purpose, since we are not attempting to derive a comprehen- 
sive isotopic composition. A similar approach was followed by 
Webber, Lee & Gupta (1992). A more complex reaction scheme 
would not in any way change our conclusions. 

We model convection in a simple way, taking a linear in- 
crease of velocity with z. Detailed self-consistent models 



of cosmic -ray driv en MHD winds ( Zirakashvili et al. 1996 



Ptuskin et al. 1997) provide explicit predictions for the convec- 



tive transport of cosmic -rays, and our approach could be used in 
future to evaluate the observational consequences of such mod- 
els. 

In this paper we concentrate on the evaluation of the B/C 
and 10 Be/ 9 Be ratios, evaluation of diffusion/convection and 
reacceleration models, and on setting limits on the halo size. 
The B/C data is used since it is the most accurately mea- 
sured ratio covering a wide energy range and having well es- 
tablished cross sections. The 10 Be/ 9 Be ratio is used rather than 
1C 'Be /(' Be + 9 Be) since it is less sensitive to solar modulation 
and to rigidity effects in the propagation. A re-evaluation of the 
halo size is desirabl e since new 10 B e/ 9 Be data are now avail- 
able from Ulysses ( |Connell 1998 ) with better statistics than 
previously. It is not the purpose of this approach to perform 
detailed source abundance calculations with a large network 
of reactions, which is still best done with the path-length dis- 
tribution approach ( DuVernois, Simpson & Thayer 1996| and 
references therein). Instead we use just the principal progeni- 
tors and weighted cro ss sections based on the obser ved cosmic- 
ray abundances (see Webber, Lee & Gupta 1992 ). Other key 
cosmic -ray ratios such as La Al/ Lt Al and sub-Fe/Fe are beyond 
the scope of this paper but will be addressed in future work. 



Also important are cosmic -ray gradients as derived from 
gamma rays; this provides a consistency check on the distri- 
bution of cosmic-ray sources, and we address this here. 

2. DESCRIPTION OF THE MODELS 

The models are three dimensional with cylindrical symme- 
try in the Galaxy, and the basic coordinates are (R,z,p), where 
R is Galactocentric radius, z is the distance from the Galac- 
tic plane, and p is the total particle momentum. The distance 
from the Sun to the Galactic centre is taken as 8.5 kpc. In the 
models the propagation region is bounded by R = Rh, z = Zh 
beyond which free escape is assumed. We take Rf, = 30 kpc. 
The range z/, = 1 - 20 kpc is considered since this is suggested 
by pr evious studies of radioactive nuclei (e.g., L ukasiak et al. 
1994a) an d the distribution of synchrotron radiation (P hillipps 
et al. 1981). For a given Zh the diffusion coefficient as a func- 
tion of momentum is determined by B/C for the case of no 
reacceleration; if reacceleration is assumed then the reaccelera- 
tion strength (related to the Alfven speed) is constrained by the 
energy-dependence of B/C. The spatial diffusion coefficient for 
the case of no reacceleration is taken as D xx = /3Do(p/ po) 5 ' be- 
low rigidity po, j3Do(p/ po) Sl above rigidity po, where the factor 
/3 (= v/c) is a natural consequence of a random-walk process. 
Since the introduction of a sharp break in D xx is an extremely 
contrived procedure which is adopted just to fit B/C at all ener- 
gies, we also consider the case 8\ = 82, i.e. no break, in order to 
investigate the possibility of reproducing the data in a phys- 
ically simpler way 1 . The convection velocity (in z-direction 
only) V(z) is assumed to increase linearly with distance from 
the plane (V > for z > 0, V < for z < 0, and dV/dz > for 
all z). This implies a constant adiabatic energy loss; the possi- 
bility of adiabatic energy gain (dV /dz < 0) is not considered. 
The linear form for V(z) is consistent with cosm ic -ray driven 
MHD wind models (e.g., Zirakashvili et al. 1996[ ). The veloc- 
ity at z = is a model parameter, but we consider here only 
V(0) = 0. 

Some stochastic reacceleration is inevitable, and it provides 
a natural mechanism to reproduce the energy dependence of the 
B/C ratio without an ad hoc form for the diffusion coe fficient 
( [Letaw, Silberberg & Tsao 19931 |Seo & Ptuskin 1994] H ein- 



bach & Simon 1995, |Simon & Heinbach 1996[ ). The spa- 
tial diffusion coefficient for the case of reacceleration assumes 
a Kolmogorov spectrum of weak MHD turbulence so D xx = 
(3Dq(p/ po) 6 with (5=1/3 for all rigidities. Simon and Heinbach 
(1995) showed that the Kolmogorov form best reproduces the 
observed B/C variation with energy. For the case of reaccelera- 
tion the momentum-space diffusion coefficient D pp is related to 
the spatial coefficient using the formula given by Seo & Ptuskin 
(1994) (their equation [9]), and Berezinskii et al. (1990) 



D PP D X 



4 P W 



38(4 - 8 2 )(4-S)w ' 



(1) 



where w characterises the level of turbulence, and is equal to 
the ratio of MHD wave energy density to magnetic field energy 
density. The main free p arameter in this relat ion is the Alfven 
speed v A ; we take w = 1 ( {Seo & Ptuskin 1994 ) but clearly only 
the quantity v\/w is relevant. 

The atomic hydrogen distribution is represented by the for- 
mula 

n HI (R 1 z) = n HI (R)e- ln2 ^ z ° ) \ (2) 
where rtfji(R) is taken from Gordon & Burton (1976) and zo 



In Paper II we considered only 8[ = and did not consider convection 
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from Cox, Kriigel & Mezger (1986) giving an exponential in- 
crease in the width of the HI layer outside the solar circle: 

/0.25kpc, fl<10kpc; 
Zo(W) -\ 0.083 e ° 11R kpc,fl> lOkpc. {i) 

The distribution of molecular hydrogen is taken from Bronfman 
et al. (1988) using CO surveys: 



n H2 (R,z) = n H2 (R)e' 



■ln2-fe/70 pc) 2 



(4) 



The adopted radial distribution of HI and H2 is shown in Fig- 
ure n 

For the ionized gas we use the two-component model of 
Cordes et al. (1991): 



»//// 



= 0.025 e 1 kpc v 20 W +0.2e 015 ^ l 2k P c ) 



cm" 3 . (5) 



The first term represents the extensive warm ionized gas and is 
similar to the distribution given by Reynolds (1989); the second 
term represents HII regions and is concentrated around R = 4 
kpc. A temperature of 10 4 K is assumed to compute Coulomb 
energy losses in ionized gas. 




R (kpc) 



FIG. 1 . — The adopted radial distribution of atomic, molecular and ionized 
hydrogen at z = 0. 

The He/H ratio of the interstellar gas is taken as 0.11 by 
number; there is some uncertainty in this quantity, but our 
value is consi stent with recent photospheric determinations 
(0.10 ±0.008: iGrevesse, Noels & Savuval 1996b- Helioseismo - 
logical methods ( |Hernandez & Christensen-Dalsgaard 1994 ) 
give a Helium abundance by mass of 0.242 corresponding to 
He/H = 0.08, but still with possible uncertainties due to the 
details of the models. Although the latter is perhaps the most 
accurate local determination, the uncertainty in extending the 
photospheric value to the interstellar medium over the whole 
Galaxy is large. Other uncertainties dominate the secondary 
production, for example the density of neutral and molecular 
hydrogen. In any case, even if He/H = 0.08 the influence of 
the uncertainty of He/H on the secondary production does not 
exceed 10%. 

The distribution of cosmic -ray sources is chosen to reproduce 
(after propagation) the cosmic -ray di stribution determined by 
analysis of EGRET gamma-ray data ( Strong & Mattox 1996 ). 
The form used is 



q(R,z) = qo 



R 

R^ 



Rq 0.2 kpc 



(6) 



where qo is a normalization constant, rj and £ are parameters; 
the /^-dependence has the same parameterization as that used 
for SNR by Case & Bhattacharya (1996), but we adopt differ- 
ent parameters in order to fit the gamma-ray gradient. We also 
compute models with the SNR distribution, to investigate the 
possibility of fitting the gradient in this case. We apply a cut- 
off in the source distribution at R = 20 kpc since it is unlikely 
that significant sources are present at such large radii. The z- 
dependence of q is nominal and reflects simply the assumed 
confinement of sources to the disk. 

We assume that the source distribution of all cosmic-ray pri- 
maries is the same. Meyer, Drury & Ellison (1997) suggest 
that part of the C and O originates in acceleration of C and O 
enriched pre-SN Wolf-Rayet wind material by supernovae, but 
the source distribution in this case would still follow that of 
SNR. 

The primary propagation is computed first giving the primary 
distribution as a function of (R,z,p); then the secondary source 
function is obtained from the gas density and cross sections, 
and finally the secondary propagation is computed. Tertiary 
reactions, such as n B — > U) B are treated as described in Ap- 
pendix A. The entire calculation is performed with momentum 
as the kinematic variable, since this greatly facilitates the inclu- 
sion of reacceleration. 

Full details of the propagation equation and numerical 
method used are given in Appendices A and B. The method 
encompasses nucleons, electrons and positrons. Energy losses 
for nucleons by ionization and Coulomb interactions are in- 
cluded following Mannheim & Schlickeiser (1994) (see Ap- 
pendix C.l). Details of the positron source function, magnetic 
fi eld and interstellar radiation field models were given in P aper 
II, and the energy loss formulae for electrons are given in the 
present paper in Appendix C.2. 

As an illustration of the calculations performed by the code, 
Figure shows the (R,z) distribution of primary 12 C and sec- 
ondary ,n B at 515 MeV/nucleon for a reacceleration model 
with Zh = 10 kpc. In practice we are only interested in the 
isotope ratios at the solar position, but it is worth noting the 
variations over the Galaxy, which are due to effect of the in- 
homogeneous distribution of sources and gas on the secondary 
production, fragmentation and energy losses. For comparison 
with gamma-ray data the full 3-D distributi on is of course im- 
portant and will be addressed in |Paper IV , but here only the 
radial cosmic-ray gradient from gamma-rays is considered. 

3. EVALUATION OF MODELS 

We consider the cases of diffusion+convection and diffu- 
sion+reacceleration, since these are the minimum combinations 
which can reproduce the key observations. In principle all three 
processes could be significant, and such a general model can be 
considered if independent astrophysical information or mod- 
els, for example for a Galactic wind (e.g., Z irakashvili et al. 
1996, p 3 tuskin et al. 1997 ), were to be used. Anticipating the re- 
sults, it can be noted at the outset that the reacceleration models 
are more satisfactory in meeting the constraints provided by the 
data, reproducing the B /C energy dependence without ad hoc 
variations in the diffusion coefficient; further it is not possible 
to find any simple version of the diffusion/convection model 
which reproduces B/C satisfactorily without additional special 
assumptions. 

In our calculations we use the B jC data summarized by Web- 
ber et al. (1996), fromHEAO-3 and Voyager 1 and 2. The spec- 
tra were modulated to 500 MV appropriate to this data using the 
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12 C (a) 

E — 515, MeV/nucleon 




10 B + 11 B 

E - 515. MeV/nucleon 



(b) 




FIG. 2._ The 3-D distribution of 12 C and 10 ' 1! 6 at 515 MeV/nucleon for reacceleration model with ih - 10 kpc, for Va = 20 km s . Parameters: see model 10500 
in Table 
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FIG. 3. — B/C ratio for diffusion/convection models without break in diffusion coefficient, for dV/dz = (solid lines), 5 (dotted lines), and 10 km s 1 kpc 1 (dashed 
lines). The cases sh own are <a\ = 1 k pr (b) zi, = 3 kpc, (c) zi, = IQkpr Solid lines- interstellar rati o shaded area: modulated to 300 - 500 MV. Dataijertical bars: 
HEAO-3, Voyager ( Webber et al. 1996 ), filled circles: Ulysses ( DuVernois, Simpson & Thayer 1996| : $ = 600, 840, 1080 MV). Parameters as in Table [J 
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FIG. 4. — B/Cratio for diffusion/convection models with break in diffusion coefficient, for dV/dz = (solid lines), 5 (dotted lines), and 10 km s~' kpc -1 (dash 
lines). The cases-shown are (a) zi, = 1 kpc, (b) zh — 5 kpc, (c) Zh — 20 kpc. Lower lines: interstellar ratio; upper lines: modulated to 500 MV. Parameters as in Table 111 
Data: as Figure Rl 



force-field approximation ( Gleeson & Axford 1968 ). We also 
sho w B/C values from Ulysses (D uVernois, Simpson & Thayer 
1996) for comparison, but since this has large modulation (600 
- 1080 MV) we do not base conclusions on th ese values. We 
use the measured 1Q Be/ 9 Be ratio from Ulysses (Connell 1998) 
and from Voyager-1,2, IMP-7/8, ISEE-3 as summarized by 
Lukasiak et al. (1994a). 



The source distribution adopted has T) = 0.5, £ = 1 .0 in eq. (|6j) 
(apart from the cases with SNR source distribution). This form 
adequately reproduces the small observed gamma-ray based 
gradient, for all Zh', a more detailed discussion is given in Sec- 
tion 4. 

3.1. Diffusion/convection models 
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FIG. 5. — ,0 Be/ 9 Be ratio for diffusion/convection models, for dV /dz = (solid lines), 5 (dotted lines), and 10 km s kpc~' (dashed lines). The cases shown are 
(a) Zh — 1 kpc, (b) zh = 5 kpc, (c) zii = 20 kpc. Data [mints from Lukasiak et al. (1994a) (Voyager-1,2: square, IMP-7/8: open circle, ISEE-3: triangle) and Connell 
(1997) (Ulysses): filled circle. Parameters as in Table fll 




FIG. 6. — Predicted 10 fie/ 9 fie ratio as functi on of la) Zj, fp. 



' dV /dz = 0, 5, 10 km s kpc~', (b) dV /dz for zh = 1-20 kpc at 525 MeV/nucleon corresponding to 
the mean interstellar value for the Ulysses data (Connell 1998); the Ulysses experimental limits are shown as horizontal dashed lines. The shaded regions show the 
parameter ranges allowed by the data. 



The main parameters are zh, Do, Si, 82 and po and dV jdz. 
We treat z/, as the main unknown quantity, and consider values 
1-20 kpc. The parameters of these models are summarized 

Table 1 For a given Zh we show B/C for a series of models 
with different dV jdz. 

Figure U shows the case of no break, S\ = 82', for each dV /dz, 
the remaining parameters Do, Si and po are adjusted to fit the 
data as well as possible. It is clear that a good fit is not pos- 
sible; the basic effect of convection is to reduce the variation 
of B/C with energy, and although this improves the fit at low 
energies the characteristic peaked shape of the measured B/C 
cannot be reproduced. Although modulation makes the com- 
parison with the low energy Voyager data somewhat uncertain, 
Figure || shows that the fit is unsatisfactory; the same is true 
even if we use a very low modulation parameter of 300 MV 
in an attempt to improve the fit. This modulation is near the 
minimum value for the en tire Voyager 17 year period (cf. the 
average value of 500 MV; Webber et al. 1996). The failure to 
obtain a good fit is an important conclusion since it shows that 
the simple inclusion of convection cannot solve the problem of 
the low-energy falloff in B/C. 

Since the inclusion of a convective term is nevertheless of in- 
terest for independent astrophysical reasons (Galactic wind) we 
can force a fit to the data by allowing a break in D xx {p), with 
Si ^ 82. Figure ^ shows cases with a break: here the parameters 
Do, Si, 82 and po are adjusted. In the absence of convection, the 
falloff in B/C at low energies requires that the diffusion coeffi- 

2 Note that the dependence of interaction rate on particle velocity itself is not sufficient to cause the full observed low-energy falloff. In leaky-box treatments the 
low-energy behaviour is modelled by adopting a constant path-length below a few GeV/nucleon, without attempting to justify this physically. A convective term is 
often invoked, but our treatment shows that this alone is not sufficient. 



cient increases rapidly below po = 3 GV (5\ ~ -0.6) reversing 
the trend from higher energies (82 ~ +0.6). Inclusion of the con- 
vective term does not reduce the size of the ad hoc break in the 
diffusion coefficient, in fact it rather exacerbates the problem 
by requiring a larger break 2 . 

Figure H shows the predicted and measured 10 Be/ 9 Be ratio; 
here we use the models with a break in D xx (p) since these do 
have the correct B/C ratio in the few 100 MeV/nucleon range 
where the Be measurements are available and are therefore ap- 
propriate for this comparison independently of the situation at 
higher energies. For our final evaluation we use 10 Be/ 9 Be data 
from Ulysses, which has the highest statistics. 

Figure || summarizes the limits on Zh and dV /dz, using the 
10 Be 1 9 Be ratio at the interstellar energy of 525 MeV/nucleon 
appropriate to the Ulysses data (Connell 1998). For Zh < 4 



kpc, the predicted ratio is always too high, even for no con- 
vection; no convection is allowed for such Zh values since this 
increases 10 Be / 9 Be still further. For Zh > 4 kpc agreement with 
10 Be/ 9 Be is possible provided < dV /dz < 7 km s" 1 kpc -1 . 
We conclude from Figure ||a that in the absence of convection 
4 kpc < Zh < 12 kpc, and if convection is allowed the lower 
limit remains but no upper limit can be set. It is interesting 
that an upper as well as a lower limit on zi, is obtained in the 
case of no convection, although w Be/ 9 Be approaches asymp- 
totically a constant value for large halo sizes and becomes in- 
sensitive to the halo dimension. From Figure |6|b, dV /dz < 7 
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FIG. 7. — B/C ratio for diffusive reacceleration models with zh = 5 kpc, Va 
= (dotted), 15 (dashed), 20 (thin solid), 30 km s (thick solid). Parameters as 
in Table bl In each case thejnterstellar ratio and the ratio modulated to 500 MV 
is shown. Data: as Figure W. 



FIG. 8. — B/C ratio for diffusive reacceleration models: zi, = 1 (dotted),^ 
(dashed), 10 (thin solid), and 20 kpc (thick solid). Parameters as in Table El 
In each case theJriterstellar ratio and the ratio modulated to 500 MV is shown. 
Data: as Figure pi 
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FIG. 9. — J.°Be/ 9 Be ratio for diffusive reacceleration models: (a) as function of energy for (from top to bottom) zi, = l t| 2, X 4. S, lXL 
as in Figure pi (b) as function of Zh at 525 MeV/nucleon corresporuiing to the mean interstellar value for the Ulysses data ( [Connell 1998 ); the Ulysses experimental 
limits are shown as horizontal dashed lines. Parameters as in Table j2j 



km s" 1 kpc" 1 and this figure places upper limits on the convec- 
tion parameter for each halo size. These limits are rather strict, 
and a finite wind velocity is only allowed in any case for zi, > 4 
kpc. Note that these results are not very sensitive to modula- 
tion since the predicted w Be/ 9 Be is fairly constant from 100 to 
1000 MeV/nucleon. 

Our results can be compared wi th those of other studies: 
zi, > 7.8 kpc (^reedman et al. 1980|), Zh < 3 kpc (Bloe men et 
al. 1993), and z h < 4 kpc ( Webber, Lee & Gupta 19920 . Most 
recently Lukasiak et al. (1994a) found 1 .9 kpc < Zh < 3.6 kpc 
(for no convection) based on Voyager Be data and using the 
Webber, Lee & Gupta (1992) models. We believe our new lim- 
its to be an improvement, first because of the improved Be data 
from Ulysses, second because of our treatment of energy losses 
(see Section fT2| ) and generally more realistic astrophysical de- 
tails in our model. The papers cited also did not consider the 
low-energy B/C data, which we have shown are in fact a prob- 
lem for diffusion/convection models. 

The cosmic-ray driven wind models of Zirakashvili et al. 
(1996) have values of dV /dz ~ 10 km s" 1 kpc" 1 , somewhat 
larger than our upper limits. Since their models are different 
from ours in many respects this is not significant, but suggests 
it would be useful to carry out calculations like those in the 
present paper for such models to provide a critical test of their 



viability. 



3.2. Diffusive reacceleration models 



The main parameters are Zh, Dq and va (po is arbitrary since 
6 is constant). Again we treat Zh as the main unknown quan- 
tity. The evaluation is simpler than for convection models since 
the number of free parameters is smaller. The parameters of 
these models are summarized in Table ||. Figure ^ illustrates 
the effect on B/C of varying va, from va = (no reacceleration) 
to va = 30 km s" 1 , for Zh = 5 kpc. This shows how the initial 
form becomes modified to produce the characteristic peaked 
shape. Reacceleration models thus lead naturally to the ob- 
served peake d form of B/C, as pointed out by several previous 
auth ors (e.g., |Letaw, Silberberg & Tsao 1993L S eo & Ptuskin 
1994, |Heinbach & Simon 1995[ )l 

Figure [| shows B/C for Zh =1—20 kpc. Our value of va ~ 20 
km s" 1 is consistent with the value obtained by Seo & Ptuskin 
(1994) which they also derived from B/C; since for stable nu- 
clei the leaky-box and diffusion treatments are equivalent this 
is a good test of the operation of our code. The value of va is 
typical of the warm ionized phase of the interstellar gas ( Seo 
& Ptuskin 1994). The exact low-energy form of B/C depends 
on details of the modulation so that an exact fit here is not 
attempted; note however that v A and Do can be (and indeed 
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Table 1 

Parameters of diffusion/convection models. 



Model 


Zh 


Do 


PO 


Si 


<5 2 


dV/dz 




kpc 


i rv28 ~ m 2 „— 1 
1U L1T1 S 


(J V 






km s ^ kpc ' 


01000 


1 


0.7 


3 


0.60 


0.60 





01010 


1 


0.7 


3 


0.60 


0.60 


10 


01 fwfl 

UlUZU 


1 


n 7 
u. / 


J 


O AO 

u.ou 


O AO 

u.ou 


90 
ZU 


03000 


3 


2.0 


3 


0.60 


0.60 





03010 


3 


1.4 


3 


0.65 


0.65 


10 


0^090 
UjUZU 


J 


i i 
i . i 


3 


O 70 
U. I\) 


O 70 
U. /U 


90 
ZU 


10000 


10 


5.0 


3 


0.60 


0.60 





10010 


10 


2.5 


3 


0.70 


0.70 


10 


1 0090 
1UUZU 


i n 
j u 


i . i 




O QO 
U.VU 


O QO 

u.yu 


90 
ZU 


01100 


1 


0.9 


5 


-0.60 


0.60 


() 


01105 


1 


0.8 


5 


-0.60 


0.60 


5 


01110 
Ul 1 1U 


1 


s 

U.o 


C 

J 


O £0 

—u.ou 


O AO 

u.ou 


1 
1 u 


03100 


3 


2.5 


5 


-0.60 


0.60 


() 


03105 


3 


2.2 


5 


-0.60 


0.60 


5 


OT1 1 O 


J 


9 n 
z.u 


C 

J 


O AO 

—u.ou 


n ao 
U.ou 


1 
1U 


04100 


4 


3.5 


5 


-0.60 


0.60 





04105 


4 


2.7 


5 


-0.60 


0.70 


5 


041 1 O 


/i 

•4 


9 S 
Z.J 


C 

J 


O AO 

—u.ou 


O 70 
U. / U 


1 
1 u 


05100 


5 


4.5 


5 


-0.60 


0.60 





05105 


5 


3.2 


5 


-0.60 


0.70 


5 


05110 


5 


2.5 


5 


-0.60 


0.70 


10 


10100 


10 


7.0 


5 


-0.60 


0.60 





10105 


10 


3.8 


5 


-0.60 


0.80 


5 


10110 


10 


3.0 


5 


-0.60 


0.80 


10 


15100 


15 


9.0 


5 


-0.60 


0.60 





15105 


15 


3.8 


5 


-0.60 


0.80 


5 


15110 


15 


3.0 


5 


-0.60 


0.80 


10 


20100 


20 


9.0 


5 


-0.60 


0.60 





20105 


20 


3.8 


5 


-0.60 


0.80 


5 


20110 


20 


3.0 


5 


-0.60 


0.80 


10 



Table 2 

Parameters of diffusive reacceleration models? 



Best fit 


Models with 


Models with SNR 


Zh 


Da 


VA 


models b 


no energy losses b 


source distribution 


kpc 


10 28 cm 2 s~ l 


km s~' 


01500 


01510 


01511 


1 


1.7 


20 


02500 


02510 


02511 


2 


3.2 


20 


03500 


03510 


03511 


3 


4.6 


20 


04500 


04510 


04511 


4 


6.0 


20 


05500 


05510 


05511 


5 


7.7 


20 


10500 


10510 


10511 


10 


12 


20 


15500 


15510 


15511 


15 


15 


20 


20500 


20510 


20511 


20 


16 


18 


Effect of varying V4 : 












05501 






5 


7.7 





05502 






5 


7.7 


5 


05503 






5 


7.7 


10 


05504 






5 


7.7 


15 


05505 






5 


7.7 


20 


05506 






5 


7.7 


25 


05507 






5 


7.7 


30 



a For all reacceleration models po = 3 GV, 8=1/3 (see Section 2 for details) 
'Parameters of the source distribution (eq. [^): r\ = 0.5, § = 1.0 
c Parameters of the SNR distribution (eq. [H]): r\ = 1.69,? = 3.33 
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were) determined from the high-energy B/C alone, and the 
low-energy agreement is then satisfactory 3 . Figure || shows 
w Be/ 9 Be for the same models, (a) as a function of energy 
for various Zh, (b) as a function of Zf, at 525 Me V/nu :leon 
corresponding to the Ulysses measurement. Comparing 
the Ulysses data point, we conclude that 4 kpc < Zh 
kpc. Again the result is not very sensitive to modulation 
the predicted 10 Be/ 9 Be is fairly constant from 100 to 
MeV/nucleon. 



with 
< 12 
since 
1000 




Kinetic energy, MeV/nucleon 

FIG. 10. — t0 Be/ 9 Be ratio for diffusive reacceleration model, showing in- 
fluence of energy losses, for zn = 1 kpc (dotted lines), 5 kpc (solid), 10 kpc 
(dashed). In each casejipper curve is with energy Jpsses, lower curve without. 
Parameters as in Table bl Data points as in Figure pi 

Figure |Io| illustrates the importance of energy losses c n the 
10 Be/ 9 Be ratio; for reacceleration cases with Zh = 1 — 2C kpc, 
we show the ratio with and without losses. Losses attenuate 
the flux of stable nuclei much more than radioactive nucle i, and 
hence lead to an increase in 10 Be / 9 Be. The effect can be s mply 
illustrated as follows. The ionization loss rate on neutral 

1.8 x 10" 7 Z 2 (n H ) f3~ l eV s" 1 , where (3 = v/c is the nubleon 



speed, and (n# ) is the average interstellar gas density. Th 
Z?<?-nuclei of 300 MeV/nucleon and for a gas disk with 
thickness and density 1 cm" 3 , (n#) = 0.2/zh cm" 3 , which 



is 



the loss time of ~ 3 x 10 years for Zh = 5 kpc. Coulomb 
on the ionized gas in the halo increase the losses furthe 
Figure |l3|); although the density is low the wide z-extent 
that the losses occur over large regions of the halo. For the 
Zh the diffusion time is w 4 x 10 8 years so the stable 9 Be 
nificantly attenuated. For the radioactive m Be (T1/2 =1.6 
years) the energy losses are negligible. Hence losses si 
cantly increase m Be/ 9 Be. As can be seen in Figure [Io|, 
ative effect is largest for large halos and becomes a domim 
feet only for Zh > 3 kpc. Although we illustrate this for the 
celeration case, the same effect applies to diffusion/convection 
models. Clearly if losses are ignored the predicted ratio 
too low and the derived value of Zh will be too small sir 
will have to be reduced to fit the observations. 

The p roton, Helium and positron spectra were presented 



the 



Paper II for the case Zh = 3 kpc using the same model as 



here, and the injection spectra were derived. The effect of 
ing the halo size is small for these spectra so we do not e 
that calculation to different Zh- 

4. COSMIC-RAY GRADIENTS 

An important constraint on any model of cosmic -ray propa 
gation is provided by gamma-ray data which give inforn ation 
on the radial distribution of cosmic rays in the Galaxy For 
a given source distribution, a large halo will give a smaller 

nova 
IWeb- 



is for 
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gives 
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(see 
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sig- 
< 10 6 
gnifi- 
rel- 
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reac- 
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be 

ce Zh 



in 
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Ktend 



ber 1997 
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theless it 
halo size 
tion. For 
Bhattach irya 
a steep ft 
Figure 
distributi 
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deduced 
MeV) 
the 7r°- 
ysis by 
gives a 
ing the 



from 
de cay i 



sizes; in 
distributi m 
lor, Man 
steeper ft 
that 
and a 
tory fit. I 
with mor 
closely tl 
ing 
not 

gradient 
beyond a 
of the 
tion reg 



prese nt 



Based 
the 

expected 
tion of 
by others 
In view 
is perhaf 
sources 
the best 



cosmic-ray gradient. It is generally believed that supe 
remnants (SNR) are the main sources of cosmic rays (see 

3 Since we are considering a ratio at the same rigidity the effect of modulation is confined to a deceleration fa 200 MeV/nucleon (cf. spectra where absolute 
intensity changes are important). 



for a recent review), but unfortunately the distribu- 
is poorly known due to selection effects. Never- 
is interesting to compare quantitatively the effects of 
on the gradient for a plausible SNR source distribu- 
illustration we use the SNR distribution from Case & 
(1996), which is peaked at R = 4-5 kpc and has 
lloff towards larger R. 

[fl] shows the effect of halo size on the resulting radial 
an of 3 Ge V cosmic-ray protons, for the reacceleration 
For comparison we show the cosmic-ray distribution 
by model-fitting to EGRET gamma-ray data (> 100 
Strong & Mattox (1996), which is dominated by 
component generated by GeV nucleons; the anal- 
Wunter et al. (1997), based on a different approach, 
sipiilar result. The predicted cosmic -ray distribution us- 
NR source function is too steep even for large halo 
fact the halo size has a relatively small effect on the 
Other related dis tributions such as pulsars ( Tay- 
hester & Lyne 1993, |Johnston \99A ) have an even 
lloff. Only for Zh = 20 kpc does the gradient approach 
obselrved, and in this case the combination of a large halo 
sli| htly less steep SNR distribution could give a satisfac- 
or diffusion/convection models the situation is similar, 
convection tending to make the gradient follow more 
e sources. A larger halo (zi, ^S> 20 kpc), apart from be- 
excluded by the l0 Be analysis presented here, would in fact 
impr )ve the situation much since Fig. |ll] shows that the 
ipproaches an asymptotic shape which hardly changes 
certain halo size. This is a consequence of the nature 
dilffusive process, which even for an unlimited propaga- 
still retains the signature of the source distribution. 



E = 3146 MeV 
z — 0.0 kpc 

model 
15-001511 
15-00251 1 
15-005511 
15-010511 
15-020511 




10 20 30 

Galactocentric radius, kpc 

FIG. 1 1 .4- Radial distribution of 3 GeV protons at z = 0, for diffusive reac- 
celeration r lodel with halo sizes zh — 1, 3, 5, 10, 15, and 20 kpc (solid curves). 
The source distribution is that for SNR given by Case & Bhattacharya (1996), 
shown as a dashed line. The msmir-ray distribution deduced from EGRET 
>100 Me\ gamma rays (Strong & Mattox 1996) is shown as the histogram. 
Parameters as in Table Q 



on these results we have to conclude, in the context of 
models, that the distribution of sources is not that 



from the (highly uncertain: see Green 1991 ) distribu- 
Slfj R. This conclusion is similar to that previously foun d 
( |Webber, Lee & Gupta 1992j |Bloemen et al. 1993] ). 
(|)f the difficulty of deriving the SNR distribution this 
s not a serious shortcoming; if SNR are indeed CR 
then it is possible that the gamma-ray analysis gives 
estimate of their Galactic distribution. Therefore in our 
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standard model we have obtained the source distribution empir- 
ically by requiring consistency with the high-energy gamma- 
ray results. 

Figure [j~2| shows the source distribution adopted in the 
present work, and the resulting 3 Ge V proton distribution, again 
compared to that deduced from gamma rays. The gradients 
are now consistent, especially considering that some system- 
atic effects, due for example unresolved gamma-ray sources, 
are present in the gamma-ray results. 

Measurements of cosmic -ray anisotropy in the 1-100 TeV 
range provide an indep endent argument for reacceleration (e.g., 
Seo & Ptuskin 1994) since the slower increase of diffusion 



coefficient with energy avoids the large anisotropics predicted 
by non-reacceleration models. Our models reproduce this be- 
haviour, the reacceleration models giving anisotropics ~ 10~ 3 
at 1 TeV, while the non-reacceleration models give > 10~ 2 . The 
observed values (~ 10~ 3 ) largely reflect the local structure of 
the interstellar magnetic field in the part of the Galaxy near 
the Sun, and hence do not give useful constraints on the large- 
scale prop agation which our model addresses (see B erezinskii 
et al. 1990). In particular it is not possible to test the large-scale 
cosmic -ray gradients at such energies by this method. It is suffi- 
cient to note that the reacceleration models are consistent with 
the observations while the non-reacceleration models are not, 
in accord with previous authors' conclusions. 



E = 3146 MeV 
z = 0.0 kpc 

model 
1 5-001500 
15-002500 
15-005500 
15-010500 
15-020500 




20 kpc 



10 20 30 

Galactocentric radius, kpc 



reacceleration is 4 kpc < zi, < 12 kpc. Our new limits should 
be an improvement on previous estimates because of the more 
accurate Be data, our treatment of energy losses, and the inclu- 
sion of more realistic astrophysical details (such as, e.g., the gas 
distribution) in our model. 

In case reacceleration is not important, the halo size limits are 
still 4 kpc < zi, < 12 kpc for the case of no convection, while 
only the lower limit holds if convection is allowed. The upper 
limit on the convection velocity gradient is dV jdz < 7 km s" 1 
kpc" 1 , and this value being allowed for large Zh only. 

The gradient of protons derived from gamma rays is smaller 
than expected for SNR sources, the closest approach to consis- 
tency being for Zh = 20 kpc; we therefore adopt a flatter source 
distribution in order to meet the gamma-ray constraints. 

The anisotropy at ~ 1 TeV predicted by our reaccelera- 
tion models is consistent with observations, while the non- 
reacceleration model predict a larger value than observed. This 
reflects the general property of such models (e.g., Seo & 
Ptuskin 1994). The large-scale propagation is however not sig- 
nificantly constrained by anisotropy measurements at the ener- 
gies considered in this paper, since local interstellar effects may 
dominate. 

The large Zh values found here would have very significant 
implications for gamma rays at high galactic latitudes, giving 
a larger inverse-Compton intensity than normally c onsidered. 



Gamma-rays will be addressed in detail in Paper IV 



We are grateful to the referee for useful suggestions. We 
thank Dr. J. J. Connell for help with the Ulysses Be data and 
for providing data prior to publication. We thank Dr. D. Bre- 
itschwerdt and Dr. V. Ptuskin for useful discussions. 



FIG. 12. — Radial distribution of 3 GeV protons at z = 0, for diffusive reac- 
celeration model with various halo sizes z.h = 1> 3, 5, 10, 15, and 20 kpc (solid 
curves). The source distribution used is shown as a dashed line, and is that 
adopted to reprodu ce the msmin-ray distrib ution deduced from EGRET >100 
MeV gamma ra^s ( Strong & Mattox 1996[ ), shown as the histogram. Parame- 
ters as in Table 



5. CONCLUSIONS 

We have shown that simple diffusion/convection models have 
difficulty in accounting for the observed form of the B/C ratio 
without special assumptions chosen to fit the data, and do not 
obviate the need for an ad hoc form for the diffusion coefficient. 
On the other hand we confirm the conclusion of other authors 
that models with reacceleration account naturally for the energy 
dependence over the whole observed range, with only two free 
parameters. Combining these results points rather strongly in 
favour of the reacceleration picture. In this case va ~ 20 km 
s , with little dependence on Zh- 

For the first time 10 Be/ 9 Be has also been computed with 
reacceleration. We take advantage of the recent Ulysses Be 
measurements to improve limits on the halo size. We empha- 
size the crucial importance of the treatment of energy losses 
in the evaluation of the 10 Be/ 9 Be ratio. The halo height with 
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APPENDIX 
A. PROPAGATION EQUATION 



The propagation equation is written in the form: 



dip x •e? /r> V7 / T7 /•> ® 2r-> ^ 1 ^ T • 



of op op p l op 



^-|(V-V)V 



T/ T> 



where = tf;(r,p,t) is the density per unit of total particle momentum, ip(p)dp = 4irp 2 f(p) in terms of phase-space density f(p), 
q(r,p) is the source term, D xx is the spatial diffusion coefficient, V is the convection velocity, reacceleration is described as diffusion 
in momentum space and is determined by the coefficient D pp , p = dp/dt is the momentum loss rate, t/ is the time scale for 
fragmentation, and r r is the time scale for the radioactive decay. The details of the numerical scheme is described in Appendix B. 

We use particle momentum as the kinematic variable since it greatly facilitates the inclusion of the diffusive reacceleration terms. 
The injection spectrum of primary nucleons is assumed to be a power law in momentum for the differ ent species, dq(p)/dp oc p~ r 
for the injected t/en^/fyn, as expected for diffusive shock acceleration (e.g., blandford & Ostriker 1980 ); the value of T can vary with 



species. The injection spectrum for 12 C and l6 was taken as dq(p)/dp oc p~ L , for the case of no reacceleration, and p~ 225 with 
reacceleration. These values are consistent with Engelmann et al. (1990) w ho give an injection index 2.23 ± 0.05. The same indices 
reproduce the observed proton and 4 He spectra, as was shown in Paper II. For primary electrons, the injection spectrum is adjusted 



to reproduce direct measurements, gamma-ray and synchrotron data; details are given in the other papers of this series (I, II, IV). 

For secondary nucleons, the source term is q(r,p) = Pcip p (r,p)[af 1 \p)nH(r) + <Ju e (p) n He(x)\ where er#(p), Oji e (p) are the pro- 
duction cross sections for the secondary from the progenitor on H and He targets, ip p is the progenitor density, and tin, nu e are the 
interstellar hydrogen and Helium number densities. 

To compute B/C and 10 Be / 9 Be it is sufficient for our purposes to treat only one principal progenitor and compute weighted cross 
sections based on the observed cosmic-ray abundances, which we took from Lukasiak et al. (1994b). Explicitly, for a principal pri- 
mary with abundance I p , we use for the production cross section a ps = J^i &' S U/Ip, where er", 7, are the cross sections and abundances 
of all species producing the given secondary. For the case of Boron, the Nitrogen progenitor is secondary but only accounts for 
10% of the total Boron production, so that the approximation of weighted cross sections is sufficient. 

For the fragmentation cross sections we use the formula given by Letaw, Silberberg & Tsao (1983). For the secondary produc - 



tion cross sections we use the Webber, Kish & Schrier (1990) and Silberberg & Ts ao (1990, see also Garcia-Munoz et al. 1987) 



parameterizations in the the form of code obtained from the Transport Collaboration (Guzik et al. 1997). Comparison of the results 
from these different versions of the cross sections gives a useful estimate of the uncertainty from this source. For the important 
B/C ratio, we take the I2 C, l6 — > 10 B, 10 C, n B, n C cross sections from the fit to experimental data given by Heinbach & Simon 
(1995). Since for Be the values of the cross sections are particularly important we give for reference the values actually used for the 
abundance-weighted cross sections at 500 MeV/nucleon, including interstellar He: a( 12 C — ► 9 Be) = 18.2 mb, a( 12 C — > 10 Be) = 8.6 
mb. For radioactive decay, r r = 771/2/I112, where t\/2 = 1.6 x 10 6 years for 10 Be. 

For electrons and positrons the same propagation equation is valid when the appropriate energy loss terms (ionization, bremsstrahlung 
inverse Compton, synchrotron) are used. Since this paper is intended to complete the description of the full model, we include the 
formulae for these los s mechan isms in Appendix C.2. A detailed description of the source function for secondary electrons and 



positrons was given in Paper II 



B. NUMERICAL SOLUTION OF PROPAGATION EQUATION 

The diffusion, reacceleration, convection and loss terms in eq. (A.l) can all be finite-differenced for each dimension (R,z,p) in the 
form 

dt~ At ~ At qi ' UJ 

where all terms are functions of (R,z,p). 



In the Crank-Nicholson implicit method (Press et al. 1992) the updating scheme is 



The tridiagonal system of equations, 



C A ' = ^ + a l ^'-a 2 ^ +A ' + a^' + q i At . (2) 
- a, V^ A ' + ( 1 + a 2 ) Vf A ' - a 3 Ci A ' = $ + 9 At , (3) 



is solved for the ip' by the standard method (Press et al. 1992). Note that for energy losses we use 'upwind' differencing to 



enhance stability, which is possible since we have only loss terms (adiabatic energy gain is not included here). 
The three spatial boundary conditions 

iP(R,z h ,p) = Tp(R-Zh,p) = mh,z,p) = (4) 

are imposed at each iteration. No boundary conditions are imposed or required at R = or in p. Grid intervals are typically AR = 1 
kpc, Az = 0.1 kpc; for p a logarithmic scale with ratio typically 1.2 is used. Although the model is symmetric around z = the 
solution is generated for — Zh < z < Zh since this is required for the tridiagonal system to be valid. 

Since we have a 3-dimensional {R,z,p) problem we use 'operator splitting' to handle the implicit solution, as follows. We apply 
the implicit updating scheme alternately for the operator in each dimension in turn, keeping the other two coordinates fixed. To 
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account for the substeps and 4- are used instead of q,, 1/r. The coefficients of the Crank-Nicholson scheme we use are given in 
Table | 

The method was found to be stable for all a, and this property can be exploited to advantage by starting with a ^> 1 (see below). 
The standard alternating direction implicit (ADI) method, in which the full operator is used to update each dimension implicitly in 
turn, is more accurate but was found to be unstable for a > 1. This is a disadvantage when treating problems with many timescales, 
but can be used to generate an accurate solution from an approximation generated by the non-ADI method. 

A check for convergence is performed by computing the timescale g ^ gt from eq. (A.l) and requiring that this be large compared 
to all diffusive and energy loss timescales. The main problem in applying the method in practice is the wide range of time-scales, 
especially for the electron case, ranging from 10 4 years for energy losses to 10 9 years for diffusion around 1 GeV in a large halo. Use 
of a time step Af appropropriate to the smallest time-scales guarantees a reliable solution, but requires a prohibitively large number 
of steps to reach the long time-scales. The following technique was found to work well: start with a large Af appropriate for the 
longest scales, and iterate until a stable solution is obtained. This solution is then accurate only for cells with a <C 1 ; for other cells 
the solution is stable but inaccurate. Then reduce Af by a factor (0.5 was adopted) and continue the solution. This process is repeated 
until a€l for all cells, when the solution is accurate everywhere. It is found that the inaccurate parts of the solution quickly decay 
as soon as the condition a < 1 is reached for a cell. As soon as all cells satisfy a < 1 the solution is continued with the ADI method 
to obtain maximum accuracy. A typical run starts with Af = 10 9 years and ends with Af = 10 4 years for nucleons and 10 2 years for 
electrons performing ~ 60 iterations per Af . In this way it is possible to obtain reliable solutions in a reasonable computer resources, 
although the CPU required is still considerable. All results are output as FITS datasets for subsequent analysis. 

More details, including the software and datasets, can be found on the World Wide Web (address available from the authors). 



B.l. DIFFUSION IN R 
As an example, the coefficients for the radial diffusion term are derived here. 

2 D xx 



RdR\ n dR I RiRi+i-Ri-i 



"'+1 T; 7T — 



Ri+i-Ri Ri-Ri-i 

Setting Rj+i = i?,— = AR, one can obtain the following expressions in terms of our standard form (eq. [A.l]) 



(5) 



ai 2Rj - AR a 2 

■=D rr .'. . ~A=D 



2Ri 



a 3 IRi+AR 
= D - 



At xt 2fl,(Afl) 2 ' Af Ri(AR) 2 ' Af AA 2fl,(A/?) 2 



(6) 



4 This corresponds to an injected flux dF(p)/dp oc /3p r or dF(E/ t )/dEi c oc p r , a form often used (e.g., Engelmann et al. 1990). Since observations are usually 
quoted as aflux with kinetic energy per nucleon as the kinematic variable a conversion is made before comparison with data: dFiE^/dE^ = -^ftip(dp j dE^) = ^-Aif) 
since dp/dE^ =A/f3, where A is the mass number, E^ is the kinetic energy per nucleon, j3 = v/c. 



Table 3 

Coefficients for the Crank-Nicholson method. 



Process 



Coordinate 



ai/Af 



a 2 /At 



a 3 /At 



Diffusion 



R 

z 



**2Ki(AR) 2 

ZV(A Z ) 2 



Uxx R,{AR? 

2ZW(Az) 2 



" 2R,(A/i) 2 

ZW(Az) 2 



Convection 11 



Z > (V > 0) 
Z < (V < 0) 
p (dV/dz > 0) 



V(Zi-i)/Az 





Vfo)/Az 
-V(z,)/Az 

-l„.4V/pi 
3 P> dz I i-l 





-V( ZM )/Az 

-In- ,£L /pi+l 



Diffusive 
reacceleration a 



2 , D P P M Dpp.i-l , 
FT /»+' i» , ' 



"PPM ( 1 2_.. 



Energy loss a 

Fragmentation 

Radioactive 
decay 



P 

R,z,p 
R,z,p 



p,/pr l 

l/3r/ 
1/3tv 








-P<=p,- Pj 
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B.2. DIFFUSIVE REACCELERATION 
In terms of 3-D momentum phase-space density f(p) the diffusive reacceleration equation is 



df(p) 
dt 



= V p .[D pp V p f(p)] = ^ 



i n 9f(p) 
P Dpp-q^- 



(7) 



The distribution is assumed isotropic so f(p) = f(p) where p = \p\. First we rewrite the equation in terms of ip(p) = 4irp 2 f(p) instead 
of f(p) and expand the inner differential: 



dtp 
~dt 



d_ 

dp 



2n d 4> 

P 'dpp~ 2 



d_ 

dp pp 



The differencing scheme is then 



D 



ppj+i 



tp M - ipt 2tp M 



Pi+i-Pi 



Pi+i 



-D 



pp,i-i 



Pi+\~Pi-\ 

In terms of our standard form (eq. [A. 1]) the coefficients for reacceleration are 

2£> WV '-i 



dip 2tp 
dp p 



ipi-ipi-i 2%pj-i 
Pi-Pi-i Pi-i 



(8) 



(9) 



O] 

At 
At 
At 



Pi+i~Pi-\ 
2 

Pi+i-Pi-i 
2Dp P j+\ 



1 2 
+ 

Pi-i Pi-i 



D p p 

Pi+i-Pi 
1 



Dpp.i-i 

Pi-Pi-I 

2 



(10) 



Pi+i-pi-i \Pi+\~Pi Pi+i 
C. ENERGY LOSSES 



For nucleon propagation in the ISM the losses are mainly due to ionization, Coulomb scattering, fragmentation, and radioactive 
decay. For electrons the important processes are ionization, Coulomb scattering, bremsstrahlung in the neutral and ionized medium, 
as well as Compton and synchrotron losses. Although all these processes are well-known the formulae for the different cases are 
rather scattered throughout the literature and hence for completeness we summarize the formulae used below. 

Figure pj| illustrates the energy loss time scales, E(dE /dt)~ x , for electrons and nucleons in pure hydrogen. The losses are shown 
for equal neutral and ionized gas number densities tin = nun = 0.01 cm" 3 , and equal energy densities of photons and the magnetic 
field U = Ub = 1 eV cm" 3 (in the Thomson limit). These gas and energy densities are chosen to characterize the average values seen 
by cosmic -rays during propagation. 




FIG. 13. — Energy loss time-scales of (a) nucleons and (b) electrons in neutral and ionized hydrogen. The curves are computed for gas densities nu = nun = 0.01 
cm -3 , and equal energy densities of photons and magnetic field U = Ub = 1 eV cm -3 (in the Thomson limit). In panel (a) solid lines correspond to ionization losses 
and dashed lines to Coulomb losses (the dotted line is that for protons). 



C.l. NUCLEON ENERGY LOSSES 



The Coulomb collisions in a completely ionized plasma a re dominated by scattering off the thermal electrons. The corresponding 
energy losses are given by (Mannheim & Schlickeiser 1994, their eqs. [4. 16], [4. 22]) 



dE\ 
dt J 



-Auric m e c 2 Z 2 n e In A ■ 



Coul 



a 2 



(1) 
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where r e is the classical electron radius, m e is the electron rest mass, c is the velocity of light, Z is the projectile nucleon charge, 
P = v/c is the nucleon speed, n e is the electron number density in plasma, x m = (3^/n/4) 1 ^ x \j2kT e /m e c 2 , and T e is the electron 
temperature. The Coulomb logarithm in the cold plasma limit is given by (e.g., Dermer 1985 ) 



InA^Ilnp^.J^ 

2 \nr e n 2 c 2 n e M+2jm e 



(2) 



where ft is the Planck constant, M is the nucleon mass, and 7 is the nucleon Lorentz factor. For the appropriate number density, 
n e ~ 10 -1 - 10~ 3 cm" 3 , and total energy E ~ 10 3 - 10 4 MeV, the typical value of the Coulomb logarithm In A lies within interval 
~ 40 - 50, instead of usually adopted value 20. 



For the ionization losses we use a general formula (Mannheim & Schlickeiser 1994, their eq. [4.24]) 

1 



dE 
~di 



(l3>l3o) = -27Tr z e cm e c 2 Z 2 

1 s=H.He 



, ^ n s [B s + B'(a f Z/(3)] 



(3) 



where a/ is the fine structure constant, n s is the number density of the corresponding species in the ISM, [3o = \Ae 2 jftc = 0.01 is the 
characteristic velocity determined by the orbital velocity of the electrons in hydrogen, and 



In 



/ 2m c c 2 /3 2 7 2 e ma 



Z s 



(4) 



where 7 is the Lorentz factor of the ion. The largest possible energy transfer from the incident particle to the atomic electron is 
defined by kinematics^] 

2m e c 2 (3 2 1 2 

Umsx l + [2 7 m e /M]' PJ 

where M m e is the nucleon mass, and I s denotes the geometric mean of all ionization and excitation potentials of the atom. 
Mannheim & Schlickeiser (1994) give the values = 19 eV and /# e = 44 eV. The shell-correction term C s /z s , the density correction 
term 8 S , and the B' correction term (for large Z or small (3) in eqs. (C.3),(C.4), can be neglected for our purposes. 
Fragmentation and radioactive decay are addressed in Appendix A. 

C.2. ELECTRON ENERGY LOSSES 



Ionization losses in the neutral hydrogen and helium are given by the Bethe-Bloch formula (Ginzburg 1979, p. 360) 

'(7-l)/3 2 £ 2 - 



dE\ 
dt ) 



j =-2nr 2 cm e c 2 — ^ Z s n 



s=H.He 



In 



2/, 2 



1 



(6) 



where Z, is the nucleus charge, n s is the gas number density, / s is the ionization potential (we use = 13.6 eV, In e = 24.6 eV, though 
the exact numbers are not very important), E is the total electron energy, 7 and (3 = v/c are the electron Lorentz factor and speed, 
correspondingly. 



The Coulomb energy losses in the fully ionized medium, in the cold plasma limit, are described by ( pinzburg 1979| , p. 361) 



dE 
dt 



= —2irr e cm e c Zn — 

Coul P 



In 



Em e c 2 \ 3 



l c 2 nZ ) 



(7) 



where Zn = n e is the electron number density. For an accurate treatment of the electron energy losses in the plasma of an arbitrary 
temperature see, e.g., Dermer & Liang (1989) and Moskalenko & Jourdain (1997). 



The energy losses due to ep-bremsstrahlung in the cold plasma are given by the expression ( von Stickforth 1961 ) 

2 7 2 / 87/3 [1-0.25(7- l) + 0.44935(7- 1) 2 -0.16577( 7 - l) 3 ] ,7 < 2; 



dE _ 

— ] =--a f r e cm e c Z n 



\I3- X [67 ln(2 7 )- 2 7 - 0.2900] , 



7>2. 



For the ee-bremsstrahlung one can obtain (Haug 1975 , Moskalenko & Jourdain 1997) 



where 



dE 
dt 



e c m(7*) = 8- 



= -^a f r 2 cm e c 2 Zn(3-f*Q cm (j*), 



(8) 



(9) 



*2 r 



An* 2 
■-=—+- 

37* 3 



,*2 



2+' 



7 



*2 



ln(p* + 7*) 



7* = V(7+D/2, 
P* 



V(7-D/2, 

5 note that there was a typing mistake in the denominator of the expression given by Mannheim & Schlickeiser (1994), which is corrected in our formula. 
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and the asterisk denotes center-of-mass variables. The total bre msstrahlung los ses in the ionized gas is the sum (dE/dt)Bi = 
{dE / dt) ep + (dE / dt) ee . A good approximation gives the expression (3inzburg 1979, p.408) 



dE\ 

— = -Aa f r 2 cm e c 2 Z(Z+\)nE 
dtj m 



ln(2 7 )-i 



(10) 



Br emsstrahlung energy losses in neutral gas can be obtained by integration over the bremsstrahlung luminosity (K och & Motz 



1959, see also paper IV| ) 



dE\ 

dt J BO 



= —c(3 n s \ dkk 



das 
~dk 



A suitable approximation (max 10% error near £ ~ 70 MeV) for eq. (C.ll) gives the combination (cf. eq. [CIO]) 

1" 

| —4afr^cm e c^E ln(27) — 

'dE' 



dt 



BO 



\ - n s M s 



n s Z s (Z s +l), 7 <;100; 

s=H,He 



7 > 800, 



(11) 



(12) 



(see 3inzburg 1979 , p. 386, 409), with a linear connection in between. Here M s is the atomic mass, and T s is the radiation length 
(T H ~ 62.8 g/cm z , T He ~ 93.1 g/cm 2 ). 



The Compton energy losses are calculated using the Klein-Nishina cross section (Jones 1965, Moskalenko & Jourdain 1997) 

^/ 7 (w)[S( 7 ,w,fcVS(7,w,r)], (13) 



dE irr 2 m e c 2 c 



dt 



2 7 2 /? 



where the background photon distribution, fy(ui), is normalized on the photon number density as n 1 = J duiup-f^uS), us is the energy 
of the background photon taken in the electron-rest-mass units, = wy(l ± (J), 



31 5 3 \ 11 3 1 1 

S(7,u;,fc)=uW k+ — + -+^r \n(2k+\) k-- + + =■ 

)[ 6 k 2k 2 ) 6 k 12(2Jt+l) 12(2fc+l) 2 



+Li 2 (-2k) 



-7<! ( fc+6 + 7 )ln(2Jfc+l)™ifc+ ' 



1 



6 4(2*+ 1) 12(2fc+ 1) ; 



■+2Li 2 (-2k) 



and Li 2 is the dilogarithm 

Li 2 (-2k)= 



-2k 



dx — ln(l — x) 



-1.6449341 +±ln 2 (2£+ l)-ln(2fc+ l)ln(2fe) + ^ i~ 2 (2k+ !)"'> > 0.2. 



The synchrotron energy losses are given by 



(14) 



(15) 



(16) 



where Ug = §- is the energy density of the random magnetic field. 
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